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Nonlinear wave dynamics in honeycomb lattices 
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We study the nonlinear dynamics of wave packets in honeycomb lattices, and show that, in 
quasi-lD configurations, the waves propagating in the lattice can be separated into left-moving and 
right-moving waves, and any wave packet composed of left (or right) movers only does not change 
its intensity structure in spite of the nonlinear evolution of its phase. We show that the propagation 
of a general wave packet can be described, within a good approximation, as a superposition of left 
and right moving self-similar (nonlinear) solutions. Finally, we find that Klein tunneling is not 
suppressed due to nonlinearity. 
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The vast interest in honeycomb lattices, which started 
more than 25 years ago in condensed matter by show- 
ing that electron waves obey the massless Dirac equa- 
tion [1, 2], has recently spread to numerous other fields. 
Examples range from electromagnetic waves in photonic 
crystals [3-7] and waveguide arrays [8-10], to cold atom 
in optical lattices [11-14] and more. However, despite 
having the same honeycomb-structured potential, there 
are also some very important differences between these 
various systems, mainly because the interactions be- 
tween the waves are different in nature. Namely, in 
graphene, the electrons have coulomb interaction and 
spin exchange, whereas EM waves can interact via nonlin- 
earity of different types (Kerr, saturable, etc. ), and cold 
atoms can be either bosons or fcrmions and display dipo- 
lar or nonlocal interactions. Naturally, it would be very 
interesting to study the effects of the different types of 
interactions on the phenomena associated with the linear 
(non-interacting) regime. For example, it was found that 
Klein tunneling in honeycomb lattices [15] is strongly 
suppressed by coulomb interaction [16]. Would inter- 
actions in other nonlinear systems also suppress Klein 
tunneling or perhaps some types of interactions preserve 
this extraordinary phenomenon ? 

Here, we study the dynamics of waves in honeycomb 
lattices, in the presence of Kerr nonlinearity, which ap- 
plies to photonic crystals, waveguide arrays and Bosc- 
Einstcin condensates (BEC). We focus on quasi-lD 
wavepackets: wave packets that are very wide in one 
direction and quite narrow in the other transverse di- 
rection. We find self-similar closed-form solutions for the 
nonlinear Dirac equation, i.e., solutions whose intensity 
structure remains unchanged during the propagation, ex- 
cept for a shift of the center. The spatial form of these 
solution can be completely arbitrary, as long as they are 



cither left movers or right movers only. Moreover, we 
show that the propagation of a general wave packet in 
the honeycomb lattice can be described using superposi- 
tion principle, to within a very good approximation, even 
in the presence of significant nonlinearity. Finally, we re- 
examine Klein tunneling in the presence of nonlinearity, 
and find that, as opposed to the electronic case, Klein 
tunneling is unaffected. 

For concreteness, we analyze here a honeycomb pho- 
tonic lattice displaying the Kerr nonlinearity (most com- 
mon optical nonlinearity). The paraxial propagation of a 
monochromatic field envelope ^ inside a photonic lattice 
exhibiting the Kerr non-linearity is described by 
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where Sn(x, y) is the modulation in the refractive index 
defining the lattice (Fig. 1 (a)), k m is the wave- number 
in the medium, no the background refractive index, and 
712 is the Kerr coefficient. The sign of n-i determines the 
type of nonlinearity, where ?i2 > corresponds to a fo- 
cusing nonlinearity (attractive interactions in the context 
of BEC). The term k m 5n/no is referred to as the opti- 
cal potential. It is convenient to transform the above 
equation to dimensionless form 
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where the coordinates are measured in units of 1 . 

Since the Kerr nonlinearity is local, each pseudo-spin 
(sub-lattice) component is expected to be affected only 
by its own intensity. We can now write the field as a two- 
component field, 'ft = (tpA iPb), where ipA,4>B are the 
amplitudes of the electric field on the two sub-lattices. 
By projecting on the Wannier states of the two lowest 
bands, it is possible to describe excitations close to the 
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FIG. 1. (a) The first Brillouin zone with the high symmetry 
points, (b) honeycomb lattice that has two sites in a unit cell. 

Dirac points by [17, 18]. 

id z ^ = n ^-Un^, n = diag(|^| 2 ,|^| 2 ), (3) 

where Ho describes the linear dynamics in the system, 
which can be elegantly (and rather accurately) expressed 
using the tight-binding approximation 

Ho = 1le{ip{k)}®a x +lm{ip{k)}®a y -V e ^{r)®l, (4) 

where tp(k) = ^j—^tj ex~p(iSj ■ k), tj's arc the hoping 
parameters, 8j are the vectors connecting the nearest 
neighbors, and V cx (r) is some additional external po- 
tential [10, 19, 20]. Considering uniaxial deformations 
t\ = t<i = t, t3 = jt, and expanding tp around one of the 
Dirac points (say K + ) (p(k = K + + p) ~ v x p x — iv y p y , 
we find that, to leading order, Eq.(3) reduces to the non- 
linear Dirac equation 



id z ^ = - 



(5) 

where v x = v/1 — (jy/2) 2 , v y = \/3j/2, and p is mea- 
sured in units of a -1 where a is the lattice constant. 

In what follows we consider a quasi- ID scenario, mean- 
ing that we consider wave packets that are very broad in 
one direction and narrow in the other. To leading order 
in the momentum, the Dirac Hamiltonian is isotropic and 
all the directions (with respect to the lattice) are equiv- 
alent (except for a numerical factor Vj). Therefore, we 
consider a wave packet moving in the x-dircction and 
hence put p y = 0. The resulting propagation equation is 

id z ^> = - [iv x d x <Z)(t x + V cx (x) <g> 1] \P - Uh$>. (6) 

The solutions of the linear equation are eigenstates of a x . 
Since the two components of the eigenvectors of a x differ 
only in their phases (±| = {a ±1^ , the nonlinear term is 
proportional to the identity operator (in the pseudo-spin 



space) . Therefore, the spinors solving the linear equation 
also solve the nonlinear equation. It is therefore sensible 
to look for solutions of the form 
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(f(x,z) ±f(x,z)). (7) 

Substituting (7) into (6) (setting V ex (x) = 0) we obtain 

8 z f(x, z) = Tv x d x f(x, z) + iU \f(x, z)\ 2 f(x, z). (8) 

This equation has general solutions of the form 

f{x, z) = g(x ± v x z) exp [iU\g(x ± v x z)\ 2 z] , (9) 

Where c/(£) is some arbitrary function. The 'up' state, 
|+), corresponds to a left moving solution g(x + v x z), 
whereas the 'down' state, |— ), corresponds to a right 
moving solution g{x — v g z). It is useful to think of the 
solutions, g(x ± v x z\ as linear combinations of plane 
waves. The right (left) moving solution is a combina- 
tion of waves the move to the right (left) and satisfy 
13 z = -v x p x (f3 z = v x p x ). 

This form of solutions has profound implications: 

1. The intensity of any wave packet is unchanged 
throughout propagation, except for some drift in 
its absolute position. 

2. The wave packet does not experience any broad- 
ening or narrowing as a consequence of the non- 
linearity, as generally happens in other nonlinear 
systems. 

3. The sign of the nonlinearity is irrelevant for the 
intensity structure. It affects only the phase of the 
wave packet. 

4. Since g(£) is arbitrary any additional noise is simply 
some other function g(^) which is also a self-similar 
solution. Hence the phenomenon of modulation in- 
stability is impossible in this system. That is, all 
wave packets are inherently stable. 

However, the most general wave packet is not com- 
posed of right (left) movers only. Rather, it is made of a 
superposition of the two, 
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where h±(x,0) = (f A (x,Q)±f B (x,0))/2. It would be 
very unusual if the dynamics of a general wave packet 
propagating in the nonlinear honeycomb lattice would be 
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the sum of the right-moving and left-moving self-similar 
solutions: 

* sup (a;, z) = h+(x, z)|+) + h-(x, z)\-) (11) 

where h±(x,z) are of the form (9). And indeed, even 
though the equation is nonlinear - where in general the 
sum of two solutions is not a solution, in this case such 
superposition turns out to be an excellent approximation. 
The reason is that the terms that 'spoil' the superposi- 
tion are products of h + and h_ (and their c.c). Since any 
physical wave packet has a finite width, 5x, after large 
enough propagation distance - the overlap between the 
h + (x — v x z) and h^(x + v x z) is negligible (since they 
move in opposite directions). Therefore, for z > fa a 
general wave packet evolves to a superposition of the left 
and right moving nonlinear solutions given in (9). We 
demonstrate the validity of this superposition principle 
for the nonlinear dynamics in honeycomb lattices by solv- 
ing Eq.(6) numerically with initial condition f(x,0) = 
Aff exp(-x 2 /a 2 ), g(x,Q) = M g x 2 cxp(-x 2 /<r 2 ), where 
AffjNg are normalization constants, set U = 0.5 and 
v x = v3/2, and compare to the superposition of the an- 
alytic solutions. We find that the agreement is excellent 
at large propagation distance (z> Sx) where the two so- 
lutions are spatially separated (Fig. 2 2 (b)), while even 
at short distances there is only very small discrepancy 
between the two ( z < a) (Fig. 2 (c) ). 

Next, we examine this unique wave dynamics when 
taking into account the fact that /3 Z is not really linear 
in p x . We do that by including higher order terms in 
the expansion of <p(k). We emphasize that <p(k) is not 
isotropic (even for a non deformed lattice, i.e., 7=1) be- 
yond the leading terms, and therefore one has to specify 
the direction with respect to the lattice. This anisotropy 
can be used to distinguish between excitations residing in 
the different Dirac cones, which suggests applications in 
controlling an electronic devices [21-23]. We proceed to 
study quasi ID dynamics in the x-direction, considering 
the next term the kinetic terms is 



f(p) = v x Px + 
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where the significance of the quadratic term can be tuned 
by controlling the deformation of the lattice. We empha- 
size that, as we approach 7 — ¥ 2, a qualitative difference 
arises between the different directions: v x — > whereas 
v y is finite and nonzero, as was demonstrated for Klein 
tunneling [10]. We find that, when the quadratic term 
is significant enough such that in the absence of nonlin- 
earity ( U = 0) the wave packet experiences noticeable 



1401 

120 

100 
80 
60 1 
40 1 
20 1 

o' 



(a) 



/ 



-200 -100 100 200 
x 

Intensity (z=140) 



(b) 



— analytic 
o numeric 



IS 



300 -200 -100 100 200 300 
X 

Intensity (z=35) 




100 200 300 
X 



FIG. 2. Simulated propagation of general wave packet (a), 
and comparison between the simulation and the superposition 
of the right and left moving solutions at z = 140 (b) and 
at 2 = 35 (c). Inset: the greatest discrepancy between the 
numeric and the analytic solutions. 



diffraction, the nonlinearity has its usual effects, that is, 
U > cause focusing of the beam, whereas U < en- 
hances the beam broadening. However, we emphasize 
that in honeycomb lattices there is a significant range 
of parameters where the quadratic terms are negligible 
and the exotic wave dynamics can be observed in exper- 
iments. 

Finally, we examine the effects of nonlinearity on one 
of the most remarkable phenomena associated with hon- 
eycomb lattices: Klein tunneling, where a wave packet 
tunnels into a potential step with probability 1 at normal 
incidence. It was found that in graphene, the coulomb 
interactions "strongly suppress Klein tunneling" [16]. 
Hence, it is interesting to study whether other types of 
interactions have such suppression effects or not. In this 
context, our type of nonlinearity represents not only EM 
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waves in photonic lattices, but also interacting BEC's. 
We solve Eq.(6) numerically in the presence of an addi- 
tional smooth step-like potential. The initial wave packet 
is composed of modes associated with the second band, 
and it is initially located in the region of the higher re- 
fractive index. In spite of the presence of nonlinearity, 
we find that the wave packet is entirely transmitted, i.e., 
it tunnels in to the region of lower refractive index with 
probability 1, meaning that Klein tunneling is not sup- 
pressed at all. This is in contrast to quasi-particles in 
graphene where Klein tunneling is strongly suppressed 
[16]. Moreover, even when we relax the Dirac approxi- 
mation, meaning taking into account 0(p 2 ) terms (Fig. 
3 (a)), the tunneling probability remains exactly 1 (Fig. 
3 (b) and Fig. 3 (c) ). We verify this result by calculating 
the projection of the left and right moving plane waves 
V± = \(p, ±\ip(z))\ 2 . The calculation reveals that right 
and left movers do not exchange population (Fig. 3 (d)). 
The fact that Klein tunneling at normal incidence is not 
effected by 0(p 2 ) terms was obtained already by [22], but 
that was for non-interacting waves only, whereas here we 
expand this finding to the nonlinear domain. However, 
since the "free" wave dynamics is strongly affected by 
the quadratic corrections, is it still surprising that Klein 
tunneling remains unaffected by it, and it is desirable to 
understand it. 

These findings can be explained in a fairly simple 
manner. The reflection amplitude is proportional to 
r oc (f3z,—p\W\/3 z ,p), where W includes the external 
potential step and the nonlinear term. At normal in- 
cidence, the states \(3 Z} p) and \(3 z ,—p) are |— ) and |+) 
respectively. Since the potential step is the same for both 
pseudo-spin components (there is no difference between 
the two sub-lattices), it is proportional to the identity 
operator. Moreover, since the incident wave packet is 
right (or left) moving, its two components are equal in 
magnitude, hence, the nonlinear term is proportional to 
the identity operator as well. Therefore, the reflection 
amplitude is proportional to the overlap of the states: 
r cx (— 1+), that completely vanishes. 

In conclusion, we studied quasi-lD waves' dynamics 
in honeycomb lattices in the presence of Kerr nonlin- 
earity. We have put a special emphasis on Bloch waves 



at the vicinity of the Dirac points, and found that the 
quasi- ID wave packet dynamics is very different from 
the 2D waves' dynamics studied previously [24, 25]. In 
fact, we have found an infinite number of non-diffracting 
self-similar solutions that are insensitive to noise, and 
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FIG. 3. (a) The propagation constant and an illustration of 
the initial momentum distribution (dashed line) (b) The in- 
tensity of the wave packet in the (x, z) plane. The dashed 
line represent the position of the potential step, (c) Illustra- 
tion of the potential step (red) and the initial (blue) and final 
(green) intensities, (d) The population of eigenstates of the 
initial (blue) and final wave packet. 



are therefore immune to modulation instability. More- 
over, we have shown that the most general wave packet 
is a superposition of left-moving and right-moving self- 
similar solutions, to a very good approximation. Finally, 
we reexamined Klein tunneling at normal incidence in the 
presence of nonlinearity, and found that, in contradistinc- 
tion to other systems exhibiting Columb interactions, the 
tunneling probability is unaffected by Kcrr-type nonlin- 
earity. 
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